	PRO GEM_RECFLUX,dir,ftime,Nproc,PS=PS
;Calculate and plot of Reconnected Flux (time integral of reconn. rate).
;Uses fix X-point and O-point.

leer_input,dir,N,nu,eta,eps,epse,dt,nstep,nstg,nstdat,nstk,nstp
x=2.*!pi*dindgen(N)/N & y=x
xpX = N/2 - 1 & ypX = N/4 - 1
xpO = N/2 - 1 & ypO = 3*N/4 - 1

aa=dblarr(N,N/Nproc) ; rectangulo de NxNproc
a=dblarr(N,N)
t=0.d0

apX=dblarr(ftime+1)
apXb=apX & apO=apX & apOb=apX &recflux=apX & tiempo=apX & recfluxb=apX

; For each frame (time)
for i=0,ftime do begin 
time = '000'+string(i,format='(i1)')
if i gt 9 then time = '00'+string(i,format='(i2)')
; Compute potencial a
for j=0,Nproc-1 do begin
file=dir+'a_'+time+'_00'+string(j,format='(i1)')
if j gt 9 then file=dir+'a_'+time+'_0'+string(j,format='(i2)')
openr,1,file,/f77_unformatted ; abre el archivo binario file
readu,1,t,aa
close,1
a(*,j*(N/Nproc):(j+1)*(N/Nproc)-1)=aa
endfor
tiempo(i) = t/eps
apX(i) = a(xpX,ypX) ; guarda el valor de a en el punto x 
apO(i) = a(xpO,ypO) ; guarda el valor de a en el punto O 
;apXb(k) = a(0,3*ypX) ; guarda el valor de a en el punto x (down)
;apOb(k) = a(0,ypO/3)   ; guarda el valor de a en el punto O (boundary)
endfor
apX = apX(0) - apX ; a(pX,t=0)=0
apO = apO(0) - apO ; a(pO,t=0)=0
recflux = (apX - apO)/eps

; Print time and reconnected flux
print,'======================================================================='
print,'Time= ',tiempo
print,' Rec. Flux = ',recflux
print,'======================================================================='
; Plot
if not keyword_set(ps) then window,0,xs=500,ys=500
if keyword_set(ps) then begin
set_plot,'PS'
device,filename=dir+'recflux_'+time+'.eps',bits_per_pixel=8,/encapsulated,xs=15,ys=12
endif
plot,tiempo,recflux,xtitle='Time',ytit='Reconnected Flux',xr=[0,max(tiempo)],yr=[min(recflux),1.1*max(recflux)],xsty=1,ysty=1,charsize=2
oplot,tiempo,recflux,psym=2
if keyword_set(ps) then begin
device,/close
set_plot,'X'
endif
end
